Mesoscopic approach to granular crystal dynamics 



Marcial Gonzalez, Jinkyu Yang, Chiara Daraio, and Michael Ortiz 
Graduate Aerospace Laboratories, California Institute of Technology, Pasadena, CA 91125, USA 

(Dated: November 17, 2011) 

We present a mesoscopic approach to granular crystal dynamics, which comprises a three- 
dimensional finite-element model and a one-dimensional regularized contact model. The approach in- 
vestigates the role of vibrational-energy trapping effects in the dynamic behavior of one-dimensional 
chains of particles in contact (i.e., granular crystals), under small to moderate impact velocities. The 
only inputs of the models are the geometry and the elastic material properties of the individual parti- 
cles that form the system. We present detailed verification results and validate the model comparing 
its predictions with experimental data. This approach provides a physically sound, first-principle 
description of dissipative losses in granular systems. 



I. INTRODUCTION 

Granular crystals, or highly-packed granular lattices, 
are strongly nonlinear systems that exhibit unique dy- 
namics. In particular, a one-dimensional chain of elastic 
spherical beads supports the formation and propagation 
of compact solitary waves as the result of Hertzian nonlin- 
ear contact interactions between the particles in the sys- 
tem [1] . The width of these waves is independent of their 

amplitude F^i — and their wave speed follows Vs oc F^^- 
This is in contrast to disordered granular media where 
the nature of the system is additionally driven by fric- 
tional and rotational dynamics. The dynamic response 
of granular crystals can be tuned by varying the particles' 
material and size, or by the application of precompression 
to the system [TH3]- Over the last decade, the response 
of these systems has drawn considerable attention and 
many potential applications have been studied, such as 
vibration mitigation [3] , shock and energy absorption [S]- 
[7], actuators [2, and sound focusing devices [S]. 

A typical experimental setup commonly used in the 
study of granular crystal dynamics consists of a monodis- 
persed one-dimensional chain of spherical beads guided 
by straight rails and impacted by a striker. Selected 
beads are instrumented with calibrated piezosensors [lOj 
that allow measurements of the forces acting inside the 
particles. These measurements confirmed the formation 
of compact solitary waves, and reported the presence of 
a characteristic amplitude decay, as the solitary waves 
propagate through the system (cf., e.g., [H]). 

A particle mechanics approach has been extensively 
used in the literature to model granular crystal dynam- 
ics. In this approach, dissipative effects are neglected and 
the interaction between beads is assumed to follow Hertz 
law [T1[3l[6HIOlfT2]. This model does indeed predict the 
formation of solitary waves of constant width that travel 
along the chain with a speed Vs oc F^^ , but it does not 
capture the experimentally-observed decay of the force. 
This discrepancy has recently motivated the inclusion of 
dissipative effects in the model, such as friction, plastic- 
ity, viscoelasticity, and viscous drag [2, 5. .in [TSHTS] . A 
fundamental challenge in modeling dissipative losses in 
granular crystals is the ability to capture quantitatively 



the amplitude decay and the wave shape variations ob- 
served experimentally, as the wave travels through the 
systems. This was done in jllj modeling dissipation in 
the form of an empirical discrete Laplacian in the veloc- 
ities, and relying on two phenomenological parameters 
obtained from best fitting experimental observations. In 
order to design engineering devices and materials that ex- 
ploit the unique wave dynamics of granular crystals, it is 
important to address the role of dissipation using a first- 
principle description that relies only on the knowledge of 
the particles' geometry and materials properties. 

This work is concerned with the formulation of first- 
principles predictive models of granular crystal dynam- 
ics. In particular, we investigate the role of vibrational- 
energy trapping effects in the overall dynamic behavior of 
one-dimensional chains under small to moderate impact 
velocities. To this end, we follow a mesoscopic approach 
comprised by two models whose only inputs are the geom- 
etry and elastic material properties of each particle in the 
system. The first model resolves the fine mesoscale struc- 
ture of dynamic collisions and explicitly accounts for the 
vibrational energy retained in each bead as the solitary 
wave propagates along the chain. For small to moderate 
impact velocities, the model is conservative and the in- 
clusion of permanent energy losses (such as plasticity and 
viscoelasticity) is not required. We achieve these prop- 
erties by abandoning the particle mechanics approach 
and adopting a three-dimensional finite-element model 
(i.e., a dynamic contact problem of three-dimensional de- 
formable elastic bodies that interact with one another 
over time is solved). 

The second model proposed in this work is a one- 
dimensional regularized contact model where the vibra- 
tional energy that remains trapped after impact is sub- 
sumed under the concept of a coefficient of restitution. 
For small to moderate impact velocities, the variation 
of this coefficient with the impact velocity is a geome- 
try and material dependent property that solely accounts 
for mesoscopic dynamic effects and that can be obtained 
from an experimental (cf., e.g., [55]) or numerical (cf., 
[TOl [T7] for elastic rods, and [H] for elastic disks) cam- 
paign of head-on collisions. This model is inherently 
energy-consistent and momentum-preserving. 

The concept of translational kinetic energy being 
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transferred to internal vibrational degrees of freedom 
during the collision of elastic particles dates back to 
Rayleigh For the case of two slowly colliding 

spheres, Rayleigh restricted attention to the fundamen- 
tal vibrational mode and estimated that the coefficient 
of restitution is a function of the impact velocity and 
does not vanish for identical spheres. However, an elas- 
tic sphere has a rich family of toroidal and spheroidal 
vibrational modes [20] and, therefore, it is expected that 
a fraction of energy is also stored in modes higher than 
the fundamental mode. Zippelius and co-workers have 
analytically shown that this is indeed the case for lower 
dimensional objects, such as rods and disks, and they 
developed microscopic models for the coefficient of resti- 
tution TBlJlSj. Unfortunately, the extension of their ap- 
proach to elastic spherical particles renders the problem 
intractable. Hence, a three-dimensional finite-element 
model is proposed in this work. 

II. MESOSCOPIC APPROACH TO GRANULAR 
CRYSTAL DYNAMICS 

Three-dimensional finite- element model. The 

first model presented in this work makes use of a three- 
dimensional finite-element mesh of the chain of beads. 
The contact constraint between beads is enforced using 
a penalty energy function pTl The trajectory of 

the elastodynamic problem naturally follows from Hamil- 
ton's principle 

= (5/[KE(q)-l/(q)-/c(q)] dt (1) 

where the generalized coordinates q{t) are the coordi- 
nates of the finite-element nodes in the deformed config- 
uration, KE(q) = iq-'^Mq is the kinetic energy, M is the 
mass matrix, ^(q) is the elastic energy, and /c(q) is the 
impenetrability or contact constraint. The impenetrabil- 
ity condition is enforced by penalizing the intersection 
between pairs of element faces (see [;21 and references 
therein). The elastic behavior of the beads is described 
by the strain-energy density of a neo-Hookean solid ex- 
tended to the compressible range [Hj. The resulting dy- 
namical system is then momentum and energy preserv- 
ing, and its trajectories are obtained by numerical time 
integration of ([T]). The multiple time scales — that coexist 
in time and space — and the complex dynamic contact of 
this three-dimensional system pose great challenges for 
numerical time integrators. We address these challenges 
with energy-stepping integrators which have proven to 
be effective and efficient in solving the dynamic behavior 
of a chain of elastic beads. These time integrators are 
energy-momentum conserving, symplectic, and conver- 
gent with automatic and asynchronous selection of the 
time step pTll^. 

Of particular interest to granular crystal dynamics 
are conservation of total energy — KE(q) -I- V{q) — and 
linear momentum. The kinetic energy and the linear 
momentum of the system can also be expressed as the 



contribution of all individual particles, i.e., KE^ and 
Jfe = Mk{q)k, where Mk and (q)^ are the mass and the 
mean velocity of particle k. Furthermore, the additive 
decomposition of the kinetic energy into a rigid-motion 
component — (KE)^ — ^Mfc|| (q)fc|p — and a vibrational 
contribution — VKEfc = KE/j, (KE)fe — will allow for in- 
vestigating the role of energy-trapping effects in granular 
crystals. 

For the purpose of validating the model, it is essential 
to identify and compute forces that resemble those mea- 
sured by embedded piezosensors in experimental setups, 
that is, following [TU], an averaged value F^^ of the con- 
tact forces acting on the particle. Any particle fc, except 
for the striker, has a contact force Fj, acting on its left 
and a contact force acting on its right; thus, the aver- 
aged force simplifies to F^^ = ^ (F\ — F'j,) • n, where fi is 
a unit vector aligned with the chain of beads. According 
to Newton's laws, the external force acting on particle k 
is 

Fr =F[,+F^,== Aj, = M,A(q), 

and, for a chain numbered from left to right, F^ ~ 
— F)._|_-^. Therefore, if the striker is referred as the first 
bead, F^^ can be readily obtained from all F°^*j,. 

One- dimensional regularized contact model. 
The second model presented in this work is motivated by 
a one-dimensional regularization of the three-dimensional 
contact problem previously presented. Then, for a one- 
dimensional chain of particles, we define the displace- 
ment Uk{t) — [{q)k{t) — (q)fc(O)] • n and the velocity 
Uk(t) = (q)fc(i) • fi, where k is the particle number- 
we assume other components of the displacement are 
zero. This intuitive dimensional reduction suggests re- 
casting the three-dimensional Lagrangian system into a 
one-dimensional mechanical system with forcing. Thus, 
displacements u{t) are given by the Lagrange-d'Alembert 
principle, i.e., 

= (5 / [(KE) (li) - V{u)] dt + jF (u, li) • Su dt (2) 

where F is a forcing term and is a one-dimensional 
regularized potential. It is worth noting that the forcing 
term results in an effective dissipation of the vibrational 
energy retained in the particles during the collision. 

We approximate the forcing term by a regularized con- 
tact model or compliant contact-force model. For direct 
central and frictionless impacts of two particles. Hunt and 
Crossley [13] proposed a compliant normal-force model of 
the form 

m{ui — U2) — — K7" — (q;k7")7 

where 7 ~ max{wi — U2,0} and j ~ ui — U2 are the 
penetration depth and speed, a is a damping factor, k 
is a spring constant, and m is the effective mass (i.e., 
= AI^^ + Af^^). This regularized contact model 
falls squarely within the dimensional reduction proposed 



3 



above, that is 



^(7) = 

Evidently, the potential energy is chosen in analogy with 
Hertz's theory, which is a good regularized model for the 
static contact problem of elastic bodies whose contact 
region remains small. Then, for heterogeneous pairs of 
linear elastic spheres, n = 3/2 and the spring constant 
K is well-known (cf., e.g., [12] )• The damping factor a is 
generally chosen to ensure that the energy dissipated dur- 
ing impact is consistent with the energy loss subsumed 
in a coefficient of restitution e. 

Many researchers have proposed approximate 
and exact [26] relationships between a and e (see, for 
example, [l^ for a detailed comparison of the predictions 
of these models with low speed impact measurements). 
In this work we adopt the exact solution proposed by 
Gonthier and co-workers f26|, i.e., the damping factor is 
determined from the implicit relation 



1 -I- avi 
1 — avie 



exp (awi(l -I- e)) 



(3) 



where e and vi are given |27j . It is worth noting that the 
subscript in Vi is not an index, it stands for the head-on 
impact velocity (i.e., the penetration speed at the start 
of the collision). This model is momentum-preserving 
and energy-consistent for all values of e. At low impact 
velocities, recent experimental data for steel [2B] suggest 
that the coefficient of restitution can be approximated 
by e = 1 — civ'^^ . Then, a is only a function of Vi, and 
the empirical coefRcients ci, C2 can be obtained from an 
experimental or numerical campaign of head-on collisions 
over a range of Vi. It bears emphasis that, in this work, 
e(vi) accounts for the vibrational energy retained in the 
beads the during collision. 

The application of Hunt-Crossley's regularized contact 
model to equations of motion ^ results in 
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7fc^^l + otkik) - 7fc+i(l + afc+i7fc+i) (4) 



where 7fc = max{ufe - Wfe+i, 0}, 7fc = Uk-iik+i: and ak = 
c^kivi.k) is given by Gonthier's energy-consistent model 
([3]). The value for Vi k is approximated by the largest 
attained relative velocity between adjacent particles k 
and fc + 1. Thus, the evolutionary equations for Vi^k and 
its time derivative ai k v[ 1^ are well-defined and given 
by: 

1. Loading/unloading conditions (7^ > 0): the irre- 
versible nature of the percussion is captured by 
Kuhn- Tucker conditions (i.e., ai,k > 0, jk—Vi^k < 0, 
and ai^kiik — Vi^k) = 0), together with a consistency 
condition if 7^ - v^k = (i.e., a^k ^ Uk - Uk+i if 
ilk - Uk+i > 0; or ai_fc = if ilk ~ Uk+i < 0). 

2. Touching/detaching conditions (7^. — 0): before 
and after the percussion of two beads, the value 



of Vi^k may experience a jump given by Vi^k = 7fe, 
at the first instant of contact (7^ > 0), or given by 
Vi^k = 0, after the collision (7^ < 0). Note that a-^^k 
then exists almost everywhere. 

A distinguishing characteristic of the proposed model is 
that damping factors ak are not assigned a priori to each 
pair-interaction but directly follow from the integration 
of equations of motion Q . This is achieved by addition- 
ally solving for the evolution in time of internal variables 
V{^k that account for the history-dependent nature of the 
dissipative process. In this respect, it is evident the simi- 
larity with models of plastic deformation — another path- 
dependent and irreversible phenomenon. 

The treatment of gravitational effects falls squarely 
within the framework considered above and it is readily 
achieved by including a forcing term in the equations 
of motion ([2]), that is 

^ = 5 j [(KE)(u)- V'(u)] At + j [F (u, u) + F9] • (5u dt 

(5) 

where = gMk and g is the gravitational constant. 
Finally, the model is integrated over time with variational 
integrators which accurately capture the energy behavior 
of forced mechanical systems [29] . 



III. VERIFICATION 

The verification of the models is twofold: (a) the as- 
sessment of the convergence and accuracy of the numer- 
ical solutions to the exact solutions of the models; and 




FIG. 1. (a) Finite-element mesh of a single spherical particle, 
(b) Coefficient of restitution e{vi) = 1 - 0.0247i;? '^^ obtained 
from best fitting (solid curve) a numerical campaign of head- 
on impacts of 9.525 mm stainless steel beads (symbols). Five 



are employed and e is computed from 



1(1 + 



where KEo and KE/ are the initial and final kinetic ener- 
gies of the center-of-mass motion obtained from the three- 
dimensional finite-element simulations. Two uniform tetra- 
hedral meshes are employed to assess convergence: (i) 53k 
nodes, 275k elements (o symbols), (ii) 410k nodes, 2.2M el- 
ements (o symbols), per spherical particle. The mesh referred 
as (i) is depicted in (a). 
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FIG. 2. (Color online) Dynamics of one bead (a,b,c) and of a one-dimensional chain (d,e,f) impacted by an identical 9.525 mm 
stainless steel bead. Predictions of the ID (3D) model are depicted by dashed (solid) curves. Energies of the 3D model are 
effectively compared by subtracting the vibrational components explicitly computed with the model (b,e). Exact conservation 
of linear momentum in both models (a,d) and total energy in the 3D model (b,e) are depicted by dash-dot curves. 



(b) the assessment of the accuracy of the one-dimensional 
regularized contact model as an approximation of the 
three-dimensional finite-element model. The former is 
formally addressed in [3T1 [551 US] j the latter is addressed 
by numerical experimentation in this work. Specifi- 
cally, we study the dynamics of one bead, and a one- 
dimensional chain, impacted by an identical 9.525 mm 
stainless steel bead with Vimp = 0.44 m/s. A uni- 
form finite-element mesh with 275k 4-node tetrahedral 
isoparametric elements and 53k nodes per spherical par- 
ticle is employed [see Fig. TJa)]. The coefficients of resti- 
tution e(wi) required in the simulations are determined 
from a numerical campaign of head-on impacts with the 
three-dimensional model. Figure [l[b) describes the pro- 
cedure followed to determine e(wi)7 

As indicated above, the regularized contact model is 
designed to exactly predict the coefficient of restitution 
of a single collision, as a consequence of solving Goth- 
ier's model ([s]). However, an accurate prediction of the 
dynamic behavior of the system is required to success- 
fully apply the model to granular dynamics. Thus, we 
compare the three-dimensional and the one-dimensional 
models in their predictions of the evolution of: (a) the 
linear momentum of each particle Jk obtained from the 
three-dimensional (respectively, one-dimensional) model 
as Jk = Mfc(q)fc • n (respectively, Jk = MkUk); (b) 
the kinetic (KE), potential (PE), and total (TE) ener- 



gies; (c) the external F™* and averaged F^^ forces. The 
vibrational energy of the three-dimensional model, i.e., 
^ VKEfc, is explicitly computed and subtracted from the 
system in order to compare the predictions with those of 
the one-dimensional model. 

Figure [2] presents the results of the verification. The 
accord between the evolution in time of the total energy 
of the one-dimensional model and the total energy of the 
three-dimensional model adjusted by the vibrational en- 
ergy retained in the system is noteworthy [Figs, and 
[SJe)]. It also bears emphasis the good agreement in the 
averaged contact force amplitude F,„ flQ] and contact du- 
ration [Figs.[2]jc) and[2]jf)]. Additionally, these numerical 
results showcase the exact conservation of linear momen- 
tum in both models [Figs.[2]^a) and[2|d)] and total energy 
in the three-dimensional model [Figs.jSJb) and[2]^e)]. 



IV. VALIDATION 

For validation purposes, we use a set of experimental 
observations for a one-dimensional vertical chain of 50 
stainless steel beads with diameter 9.525 mm, impacted 
by an identical stainless steel bead with several Uimp- 
The material properties of the beads are E — 200 GPa, 
v = 0.3, and the density is p = 7900 kg/m'^. The time 
history of forces measured by calibrated piezosensors in- 
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FIG. 3. (Color online) One-dimensional regularized contact 
model for Wimp = 0.63 m/s. (a) Time history of measured 
forces (dashed curves) and numerical predictions with and 
without gravitational forces, (b) Conservation of total linear 
momentum and time history of individual Jt- (c) Time his- 
tory of total (TE), kinetic (KE), and potential (PE) energies. 



serted in even-numbered beads is presented in Fig.jsj^a) — 
the bead impacted by the striker is referred as the first 
bead. Characteristically, the measured force decays as 
the solitary wave propagates along the chain. 

Figure ^a.) shows the time histories of the forces pre- 
dicted by the one-dimensional regularized contact model 
with and without gravitational effects, i.e., by equations 
of motion ^ and respectively. The good agree- 
ment with the experimental observations is evident in 
the figure. The total linear momentum of the system 
without gravitational forces is a constant of motion, and 
this momentum-preserving property of the model is ver- 



FIG. 4. (Color online) Comparison of predictions obtained 
with our one-dimensional numerical model, with and without 
gravitational forces, with measured values (empty dots) for 
all five «imp. (a) Evolution of the solitary wave full width 
at half maximum as a function of the particles position in 
the chain, (b) Decay of the maximum force, (c) Speed of 
the solitary wave versus the maximum force and its best fit 



Five different Himp are used, 0.31, 0.44, 0.63, 0.1 



and 1.25 m/s. Error bars represent the standard deviation 
from five samples of each data point. 



ified in Fig. |3[b). The decaying behavior in the total 
energy of the system — that resembles the vibrational en- 
ergy trapped in the beads — is observed in Fig. ^c). It 
is worth noting that both models presented in this work 
are intertwined in the validation since e(vi) is obtained 
from three-dimensional finite-element simulations. 

It is also interesting to note that the mesoscopic ap- 
proach predicts all the well-known qualitative behav- 
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ior of one-dimensional granular crystal dynamics, with 
and without the addition of the nonuniform gravitational 
preload. Particularly, we observe in Fig. [4j (a) the for- 
mation and propagation of solitary waves with a finite 
width (measured as full width at half maximum) that 
follows the experimental trend and that is independent 
of F^, Fig. Qa)- (b) a decay of the force that follows 
cx Fig.gb); (c) a solitary wave speed that fol- 

lows the experimentally-observed trend and Vs oc F^^^^ 
Fig. ^c). Evidently, gravitational effects are not neg- 
ligible for the length of the chain and for the impact 
velocities considered in this work. Most importantly, 
vibrational-energy trapping effects appear to play a rel- 
evant role in the dynamic behavior of one-dimensional 
granular crystals under small to moderate impact veloci- 
ties, and they explain the experimentally-observed decay 
behavior of Fm. 



V. SUMMARY AND CONCLUSIONS 

We have presented a mesoscopic approach to granular 
crystal dynamics, which comprises a three-dimensional 
finite-element model and a one-dimensional regularized 
contact model. The approach aims at investigating the 
role of vibrational-energy trapping effects in the overall 
dynamic behavior of one-dimensional chains under small 
to moderate impact velocities. The three-dimensional 
finite-element model resolves the fine mesoscale struc- 
ture of dynamic collisions and explicitly accounts for 
the vibrational kinetic energy retained in each bead as 
the solitary wave propagates along the chain. The one- 
dimensional regularized contact model accounts for meso- 
scopic dynamic effects by means of a restitution coeffi- 
cient, i.e., the vibrational energy that remains trapped af- 
ter impact is subsumed under the concept of a coefficient 
of restitution. We have specifically adopted the compli- 
ant normal-force model proposed by Hunt and Crossley, 
and the damping-factor model proposed by Gonthier and 
co-workers. The only inputs of these models are the ge- 
ometry and the elastic material properties of the indi- 



vidual particles that form the granular crystal (e.g., we 
have obtained the variation of the coefficient of restitu- 
tion with the impact velocity from a numerical campaign 
of head-on collisions) . 

We have also presented a detailed verification and val- 
idation of the mesoscopic approach, which include: (a) 
the assessment of the accuracy of the one-dimensional 
regularized contact model as an approximation of the 
three-dimensional finite-element model; (b) the one-to- 
one comparison of the experimental and simulated time 
histories of averaged forces in a one-dimensional chain of 
50 stainless steel beads impacted at five different veloci- 
ties covering the range from 0.31 to 1.25 m/s. The good 
agreement of the later and the ability of the model to 
predict well-known properties of one-dimensional granu- 
lar crystal dynamics (e.g., the formation of solitary waves 
with a finite width that is independent of the solitary 
wave amplitude, the decay of the force as the solitary 
wave propagates along the chain, and the scaling of the 
solitary wave speed with the wave amplitude) are note- 
worthy. 

Finally, we remark that the mesoscopic approach pre- 
dicts the experimentally-observed decay of the solitary 
wave amplitude from first-principles for the first time by 
solely addressing the role of vibrational-energy trapping 
effects. We thus conclude that these effects play a central 
role in the dynamic behavior of one-dimensional granular 
crystals under small to moderate impact velocities. The 
extension of this approach to multi-dimensional granular 
crystals is a worthwhile direction of future research. 
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